My Title

library(aqp)

## This is aqp 1.33

library(soilDB)
library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

library(mapview)

## Warning: package 'mapview' was built under R version 4.1.2

## Warning: multiple methods tables found for 'crop'

## Warning: multiple methods tables found for 'extend'

# get pedons from the selected set
pedons <- fetchNASIS(from = 'pedons')

## Loading required namespace: odbc

## multiple horizontal datums present, consider using WGS84 coordinates (x_std, y_std)

## Warning: 'the_value' has been rounded to the nearest integer.

## mixing dry colors ... [125 of 5482 horizons]

## mixing moist colors ... [148 of 7224 horizons]

## NOTE: some records are missing rock fragment volume

## -> QC: some fragsize_h values == 76mm, may be mis-classified as cobbles [1714 / 12902 records]

## NOTE: all records are missing artifact volume

## replacing missing lower horizon depths with top depth + 1cm ... [113 horizons]

## top/bottom depths equal, adding 1cm to bottom depth ... [15 horizons]

## -> QC: sites without pedons: 
##  Use `get('sites.missing.pedons', envir=soilDB.env) for site record IDs (siteiid)`

## -> QC: duplicate pedons: 
##  Use `get('dup.pedon.ids', envir=soilDB.env) for pedon record IDs (peiid)`

## -> QC: horizon errors detected:
##  Use `get('bad.pedon.ids', envir=soilDB.env)` for pedon record IDs (peiid)
##  Use `get('bad.horizons', envir=soilDB.env)` for horizon designations

## -> QC: pedons missing bottom hz depths:
##  Use `get('missing.bottom.depths', envir=soilDB.env)` for pedon record IDs (peiid)

## -> QC: equal hz top and bottom depths:
##  Use `get('top.bottom.equal', envir=soilDB.env)` for pedon record IDs (peiid)

# REMOVE PEDONS WITH NA COORDINATES
pedons.sp <- subset(pedons, !is.na(x_std) & !is.na(y_std))

# create sf object (more on this in next section)
pedon.locations <- st_as_sf(
  site(pedons.sp), 
  coords = c('x_std','y_std'),
  crs = st_crs(4326) #WGS84 GCS
)

# plot an interactive map
mapview(pedon.locations, 
        legend = FALSE, 
        map.types = 'OpenStreetMap', 
        label = pedon.locations$site_id)